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ABSTRACT 


The  major  thrust  of  the  research  undertaken  was  to  find 
the  most  efficient  methods  for  solving  viscous  flow  problems 
governed  by  the  N a v i e r - S t ok e s  equations.  Both  external  flows 
past  obstacles  and  internal  flows  in  tanks  and  channels  were 
investigated.  Two  aspects  in  particular  have  been  studied 
most  thoroughly.  These  are  the  f i n i t e - e 1  erne n t  discretization 
of  the  flow  equations  and  the  implementation  aspects  involv¬ 
ing  most  effective  use  of  hardware  in  order  to  obtain  small 
solution  times. 

The  first  topic  involved  a  major  mathematical  investiga¬ 
tion,  and  the  problems  have  been  solved  essentially  in  their 
entirety.  Specifically,  simple  rigorous  tests  have  been 

developed  for  testing  the  stability  of  any  potential  finite 
element  scheme.  This  problem  which  has  hitherto  resisted 
numerous  attempts  at  a  general  solution  is  now  reduced  to 
a  simple  calculation  with  a  small  number  of  variables  which 
is  carried  out  a  priori  without  the  need  for  extensive  compu¬ 
ter  calculations. 

The  second  topic,  which  reduces  to  finding  effective 
algorithms  for  solving  the  huge  systems  of  nonlinear  alge¬ 
braic  equations  yielded  by  the  discretization  methods  has 
been  subdivided  into  two  parts,  associated  with  low  and  high 
Reynolds  number  flows  respectively.  For  both  of  these  cases 
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now  methods  have  been  proposed  and  shown  to  be  effective  both 
theoretically  and  computationally. 

Effort  was  also  devoted  toa  number  of  other  problems  of 

interest.  First,  the  question  of  periodic  acoustic  wave  pro¬ 

pagation,  solved  by  mixed  finite  element  methods  was  investi¬ 
gated.  An  extension  of  the  theory  for  viscous  flows  was 

required  and  has  been  obtained.  Second,  the  stability  of 
Gaussian  elimination  applied  to  discrete  Helmholtz  problems, 
without  pivoting,  was  studied  via  model  problems.  It  was 

shown  that  contrary  to  various  conjectures  made  in  the  litera¬ 
ture  the  technique  is  unstable  even  under  'stabilizing* 

perturbations  of  the  wave  number.  Third,  it  was  realized 
that  close  connections  exist  between  the  finite  element 

theory  for  low  Reynolds  number  and  that  for  the  problems  of 

linear  e  1  a s t ome c h a n i c s  posed  in  terms  of  stresses.  The 

results  obtained  for  the  fluid  case  were  then  applied  to 
the  solution  of  some  outstanding  problems  of  the  elasticity 


1.  Na v i e r- S t ok c s  flows 

1.1  New  finite  element  classes.  [2]  [4]  [5]  [6]  [13] 

Several  new  classes  of  finite  element  pairs  have  been 

introduced,  which  are  of  optimal  accuracy  in  the  H^(-)  x 

2 

L  (...)  setting.  Some  of  these  pairs  have  been  implemented 
in  the  Stokes,  Oseen  and  N a v i e r - S t o k e s  cases.  The  results 
of  the  computer  programs  are  as  predicted  and  show  in  some 
cases  an  order  of  magnitude  improvement  in  accuracy  for  the 
same  amount  of  work  as  with  conventional  schemes. 


1.2  Stability  of  Navier-Stokes  discretization  [ 7  ]  ,  [9],  [ 1 0  ]  ,  [14] 

The  investigations  revealed  a  remarkable  dichotomy  of 
behaviour  between  low  order  and  high  order  finite  element 
approximations.  When  the  investigation  began,  the  subject 
consisted  of  a  number  of  special,  often  difficult  proofs  that 
a  small  number  of  discretizations  indeed  yielded  stable  appro¬ 
ximations  in  the  sense  that  as  some  mesh  parameter  h  tended 
to  zero,  the  computed  solutions  approached  the  exact  solution 
of  the  Navier-Stokes  equations.  The  finite  elements  covered 
by  these  proofs  by  no  means  included  those  used  in  practical 
engineering  computations.  The  first  problem  solved  consists 
of  giving  a  test  whereby  any  finite  element  pair  of  higher 
than  quadratic  order  in  velocities  can  be  relatively  simply 
processed  and  its  stability  or  otherwise  verified.  The  test 
was  applied  to  verify  the  stability  of  numerous  existing 
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computations  and  then  to  devise  better  and  stable  higher 
order  schemes.  The  existence  of  this  test  completely  dis¬ 
poses  of  the  stability  question  for  higher  order  element 
pairs.  Thus,  an  outstanding  problem  was  resolved. 

Many  engineering  calculations  in  the  aerospace  and  other 
fields  are  performed  using  element  pairs  of  lower  order  than 
quadratic.  One  such  pair,  previously  a  candidate  for  the 
'work-horse'  of  viscous  flow  calculations,  is  the  bilinear/ 
constant  combination.  Many  calculations  have  been  carried 

out  with  these  elements.  A  major  result  of  our  work  on  the 
lower  order  elements  has  been  that  the  pressure  field  for 
this  element  does  not_  converge  to  the  true  pressure  field 
of  the  N a v i e r- S t oke s  equations  as  the  meshsize  h  approaches 
zero.  Although  we  showed  how  the  resulting  field  could  be 
reprocessed  to  obtain  a  useful  pressure  field  in  simple 
cases,  it  is  not  at  all  clear  that  this  can  be  done  in 

general.  Thus,  these  popular  elements  have  a  basic  flaw. 

To  circumvent  this  difficulty,  several  element  pairs  of 
comparable  simplicity  have  been  introduced  and  proved  to  be 
stable  and  hence  useful  for  reliable  computation.  Neverthe¬ 
less,  the  fact  remains  that  low  order  elements  are  to  be 

suspected  of  inherent  instability  from  now  on,  in  the  absence 
of  rigorous  proofs  to  the  contrary. 

The  work  on  low  order  elements  was  thought  to  be  of 

sufficient  interest  to  warrant  an  abstract  theorv.  since 


similar  problems 


arise 


in  other  fields  of  computational 


mechanics.  This  work  continues. 

1.3  Solution  algorithms.  [7],  [11],  [12],  [13] 

For  the  low  Reynolds  number  case,  a  method  called  the 
"semi-direct"  solution  algorithm  has  been  introduced.  To 
explain  it,  consider  the  system  of  equations 

Au  +Bp=f  A  =  A  (  u  ) 

T 

B  u  =  g 

which  is  typical  of  the  form  assumed  by  the  stationary  Navier- 
Stokes  after  discretization  by  finite  element  methods.  Add¬ 
ing  eP ( le[  <  <  1)  to  both  sides  of  the  second  equation  and 

iterating  according  to 


Au 


n  +  1 


B  p 


f 


_T  n  +  1 
B  u 


"  e  P 


n  +  1 


n 

8  -  ep 


leads  to  the  inner  iteration 


(A  +  EBBT)  un+1  =  f  -^B(g  -  epn). 

original 
from  the 
required 


The  tremendous  advantage  of  this  problem  over  the 
one  is  that  the  pressure  pn+^  has  been  eliminated 
calculation.  Because  of  this,  computer  resources 


are  reduced  by  a  factor  of  10  or  more  over  the  original  sys¬ 


tem.  The  inner  iteration  is  implemented  by  a  fast  conjugate 
qradient  algorithm.  This  technique  has  been  thoroughly 

analysed  and  used  by  us  in  many  different  settings,  always 
performing  well,  provided  Re  is  not  too  large.  It  has  in 
fact  become  a  standard  technique  for  these  flows,  offering 
the  advantage  over  the  more  usual  penalty  method  of  comput¬ 
ing  the  exact  solution  of  the  system,  rather  than  the  exact 
solution  of  a  perturbed  system. 

For  the  higher  Reynolds  numbers,  the  semi-direct  method 
is  not  needed,  because  the  calculation  can  now  be  treated 
explicitly  (convection  dominated  rather  than  diffusion 
dominated).  Here,  another  problem,  also  associated  with  the 
pressure  is  encountered.  The  problem  is  that  of  enforcing 
the  constraint  of  incompressibility  which  takes  the  form  of 
the  second  of  the  above  equations,  namely  B  u  =  g  =  0.  A 
new  method  for  handling  this  has  been  introduced,  based  on 
the  observation  that  for  finite  element  discretizations  the 


transpose  of 

the  matrix  B 

occurs 

in 

the 

momentum  equation. 

The  technique 

is  to  observe 

that 

the 

ve 

locity  computation 

requires  the 

value  of  Bp 

only, 

and 

no 

p  itself.  But 

Bp  is  simply  the  least  squares  solution  of  B  times  the 
momentum  equation.  This  least  squares  computation  can  be 
done  very  cheaply  (one  or  two  iterations)  in  a  time  marching 
regime  and  thus  the  pressure  difficulty  may  be  overcome. 


2.1  Periodic  Acoustic  Wave  Propagation.  [1],  [3] 

An  extension  of  the  analysis  applicable  to  the  viscous 
flow  problems  is  required  in  order  for  the  results  to  be 
valid  for  the  acoustics  case.  This  analysis  has  been  success¬ 
fully  carried  out,  with  the  result  that  a  significant  improve¬ 
ment  in  accuracy  for  essentially  the  same  amount  of  work  over 
the  least  squares  method  previously  used,  is  obtained. 

A  related  study  considered  the  solution  of  periodic 
acoustic  phenomena  governed  by  the  Helmholtz  equation.  The 
objective  here  was  to  prove  or  disprove  the  proposition  that 
Gaussian  elimination  could  be  applied  directly  without 
partial  pivoting,  provided  the  wave  number  was  not  an  eigen¬ 
value  of  the  Laplacian  operator.  If  true,  this  would  consi¬ 
derably  lessen  the  computer  resources  required  for  solution. 
Unfortunately,  as  we  showed,  the  proposition  is  false.  A 
rigorous  analysis  showed  that  even  after  making  perturbations 
of  the  meshsize  order  to  move  the  eigenvalues  away  from 
singular  points  the  use  of  pivoting  is  essential  in  general. 
Otherwise,  unbounded  growth  of  the  elimination  multipliers 
can  occur,  and  the  ensuing  growth  of  the  roundoff  errors  will 
render  the  computations  meaningless. 
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